#ifndef Algebra_Space_h
#define Algebra_Space_h

#include "AlgebraModule.h"
#include "IndexArray.h"
#include "IndexMap.h"

// #include <Eigen/Dense>
// #include <Eigen/Sparse>

SG_NAMESPACE_OPEN
namespace Algebra
{
    // Forward declarations.
    class Sparsity;
    template <typename TDataType> class Space;

    /** an index set containing a list of integers. */
    typedef std::vector<Index_t> GeneralIndexSet;

    /** an index set containing a list of evenly spaced integers */
    struct StrideIndexSet
    {
        StrideIndexSet()
        {
            n = 0;
            first = 0;
            step   = 1;
        }

        /** A list of index, \f$ first : (n - 1) * step + first : step \f$ */
        StrideIndexSet (const Index_t& n_, const Index_t& first_, const Index_t& step_)
        {
            n = n_;
            first = first_;
            step  = step_;
        }

        Index_t n;      //< the length of the locally owned portion of the index set
        Index_t first;  //< the first element of the locally owned portion of the index set
        Index_t step;   //< the change to the next index
    };

    /**
     * \addtogroup LinearSpace Linear Space
     *
     * One of the main problems to use a linear solver or incorporate a new one is the matrix and
     * vector types they use which can be incompatible with those used by the application. In order
     * to solve this problem another layer of abstraction is necessary in order to encapsulate matrices
     * and vectors operations and to be used as an interface to new matrix and vector types.
     *
     * Space deﬁnes a matrix and a vector and also their operators. Encapsulating deﬁnitions and
     * operators in one object lets diﬀerent developers to implement their interfaces without any
     * conﬂict. Space deﬁnes a simple interface with all operations which are required by linear solvers.
     *
     * @see Dadvand P. A framework for developing finite element codes for multi-disciplinary applications.
     *
     * GoF's Adapter Pattern is recommented to adapt the udnerlying algebraic types to the desired interfaces.
     *
     * Convert the interface of a class into another interface clients expect. Adapter lets classes work
     * together that couldn't otherwise because of incompatible interfaces.
     *
     * Target defines the domain-specific interface that Client uses. Adaptee defines an existing interface
     * that needs adapting. Adapter adapts the interface of Adaptee to the Target interface.
     *
     * We can do this in one of two ways: (1) by inheriting Target's interface and Adaptee's implementation
     * or (2) by composing a Adaptee instance within a Adapter and implementing Adapter in terms of Adaptee's
     * interface. These two approaches correspond to the class and object versions of the Adapter pattern.
     *
     * You need to use several existing subclasses, but it's impractical to adapt their interface by subclassing
     * every one. An object adapter can adapt the interface of its parent class.
     *
     * Consider the following codes that built on the algebraic objects,
     * 
     * \code{.cpp}
     * void NewmarkBeta::Step (Vec_t &U, Vec_t &V, Vec_t &A)
     * {
     *   SpMat_t M, C, K;
     *   mODEs->ComputeImplicitJacobian (t, U_k, V_k, A_k, K, C, M);
     *   // ...
     * }
     * \endcode
     *
     * it's usually necessary to instantiate
     * the cooresponding classes inside the target' client. Put it another way, a methology shoud be introduced
     * to create the needed algebraic objects. For this, the Abstract Factory Pattern is taken to provide an
     * interface for creating families of related or dependent objects without specifying their concrete classes.
     *
     * @see Erich Gamma. Design Patterns:elements of reusable object-oriented software. Addison Wesley, 1994.
     */
    //@{

    /**
     * @brief Define the interfaces of the sequential/parallel vector.
     *
     * Vectors are used to store discrete PDE solutions, right-hand sides for linear systems, etc.
     *
     * @author nene
     * @date ‌Januray‌, 2025.
     */
    template <typename TDataType>
    class ALGEBRA_EXPORT Vector
    {
    public:
        using DataType = TDataType;

        Vector () {};
        explicit Vector (Sparsity* sp) : mSparsity(sp) {};
        virtual ~Vector () {};

        /** The corrresponding linear space which this algebraic object pertained to.*/
        virtual Space<TDataType>* getSpace () = 0;

        // /** Test whether if the vector has been initialized. */
        // virtual bool initialized () const
        // {
        //     return mInitialized;
        // }

        // /** Initialize this vector using the sparsity structure. */
        // virtual void init (Sparsity* sp) = 0;

        /**
         * Calls the Vector's internal assembly routines, ensuring
         * that the values are consistent across processors.
         */
        virtual void close () = 0;

        /**
         * Restores the \p Vector<T> to a pristine state.
         */
        virtual void clear () = 0;

        /** Start assembling the elements stashed. */
        virtual void beginAssembly () = 0;

        /** Finish assembling the elements stashed. */
        virtual void endAssembly () = 0;

        /** Perform the assembly operation. */
        virtual void assemble () = 0;

        // /** Test whether the vector has been assembled. */
        // virtual bool assembled () = 0;

        /** Resize *this to a vector of length size. */
        virtual void resize (Index_t size) = 0;

        /** Resize the vector to size while retaining old values. */
        virtual void conservativeResize (Index_t size) = 0;

        /** The number of rows of the vector */
        virtual Index_t rows () const = 0;

        /** The number of rows owned by this process. */
        virtual Index_t local_rows () const
        {
            return row_end () - row_start ();
        }

        /** The index of the first vector element stored on this processor.*/
        virtual Index_t row_start () const = 0;

        /** The index of the last vector element (+1) stored on this processor. */
        virtual Index_t row_end () const = 0;

        /** Set all entries to 0. */
        virtual void setZero () = 0;

        /** Set v(i) = \p value. */
        virtual void set (const Index_t& i, const TDataType& value) = 0;

        /** Add \p value to the vector entry specified by \p i. */
        virtual void add (const Index_t& i, const TDataType& value) = 0;

        // /** The reference of the ith entry of the vector. */
        // virtual TDataType& operator() (const Index_t& i) = 0;

        /** A copy of the ith entry of the vector.*/
        virtual TDataType operator() (const Index_t& i) const = 0;

        // /** The reference of the ith entry of the vector. */
        // virtual TDataType& operator[] (const Index_t& i) = 0;

        /** A copy of the ith entry of the vector.*/
        virtual TDataType operator[] (const Index_t& i) const = 0;

        /** Print the elements of a vector to the given output stream. */
        virtual void print (std::ostream& out) const = 0;

        /** Create a subvector using the indices in \p rows. */
        virtual void create_subvector (Vector<TDataType>& subvector, const GeneralIndexSet& rows) const = 0;

        /** Create a subvector from the stride index set in \p rows. */
        virtual void create_subvector (Vector<TDataType>& subvector, const StrideIndexSet& rows) const = 0;

        /** Determine the vector component with maximum real part \p value and its location \p index. */
        virtual void max (Index_t& index, TDataType& value) const = 0;

        /**
         * The \f$ \ell_{\infty} \f$-norm of the vector, i.e. the maximum
         * absolute value of the entries of the vector.
         */
        virtual void linfty_norm (TDataType& value) const = 0;

      protected:
        // bool      mInitialized = false;    //< whether or not the vector has been initialized
        Sparsity* mSparsity    = nullptr;  //< the data layout across processes.
    };

    /**
     * @brief Define the interface of the dense(full) matrix.
     *
     * The dense ("full") matrix, that is, a matrix where we store every element whether it is
     * zero or nonzero. It's useful for storing element stiffness matrices before summation
     * into a global sequential/parallel matrix.
     *
     * @author nene
     * @date ‌Januray‌, 2025.
     */
    // template <typename TDataType>
    // using DenseMatrix = Eigen::Matrix<TDataType, Eigen::Dynamic, Eigen::Dynamic>;

    template <typename TDataType>
    class ALGEBRA_EXPORT DenseMatrix
    {
    public:
        typedef TDataType DataType;

        /** Print the elements of a matrix to the given output stream. */
        virtual void print (std::ostream& out) const = 0;
    };

    /** A small structure to hold a non zero as a triplet (i,j,value). */
    template <typename TDataType>
    class Triplet
    {
    public:
        Triplet () : m_row (0), m_col (0), m_value (0)
        {}

        Triplet (const Index_t& i, const Index_t& j, const TDataType& v = TDataType (0)) : m_row (i), m_col (j), m_value (v)
        {}

        /** The row index of the element */
        const Index_t& row () const
        {
            return m_row;
        }

        /** The column index of the element */
        const Index_t& col () const
        {
            return m_col;
        }

        /** The value of the element */
        const TDataType& value () const
        {
            return m_value;
        }

    protected:
        Index_t      m_row, m_col;
        TDataType    m_value;
    };

    // typedef Eigen::Triplet<Real_t>    Triplet_t;
    // typedef Eigen::Triplet<Complex_t> ComplexTriplet_t;

    /**
     * @brief Define the interfaces of the sequential/parallel sparse matrix.
     *
     * sparse matrix contains contribution of local element matrices. Localization of local element
     * matrix into global (structural) matrix is determined using element code numbers.
     *
     * Using a common base class allows for uniform access to sparse matrices from various different
     * solver packages in different formats.
     *
     * @author nene
     * @date ‌Januray‌, 2025.
     */
    template <typename TDataType>
    class ALGEBRA_EXPORT SparseMatrix
    {
    public:
        typedef TDataType                          DataType;

        SparseMatrix () {};
        explicit SparseMatrix (Sparsity* sp) : mSparsity(sp) {};
        virtual ~SparseMatrix () {};

        /** The corrresponding linear space which this algebraic object pertained to.*/
        virtual Space<TDataType>* getSpace () = 0;

        // /** Test whether if the matrix has been initialized. */
        // virtual bool initialized () const
        // {
        //     return mInitialized;
        // }

        // /** Initialize this matrix using the sparsity structure. */
        // virtual void init (Sparsity* sp) = 0;

        /** Restore to a pristine state. */
        virtual void clear () = 0;

        /**
         * Call the SparseMatrix's internal assembly routines, ensuring
         * that the values are consistent across processors.
         */
        virtual void close () = 0;

        /** Start assembling the elements stashed. */
        virtual void beginAssembly () = 0;

        /** Finish assembling the elements stashed. */
        virtual void endAssembly () = 0;

        /** Perform the assembly operation. */
        virtual void assemble () = 0;

        /** Test whether the matrix has been assembled. */
        virtual bool assembled () const = 0;

        /** Resize *this to a vector of length size. */
        virtual void resize (Index_t rows, Index_t cols) = 0;

        /** Resize the vector to size while retaining old values. */
        virtual void conservativeResize (Index_t rows, Index_t cols) = 0;

        /** The global number of rows of the matrix */
        virtual Index_t rows () const = 0;

        /** The global number of columns of the matrix */
        virtual Index_t cols () const = 0;

        /** The local number of rows owned by this process. */
        virtual Index_t local_rows () const
        {
            return row_end () - row_start ();
        }

        /** The local number of columns owned by this process. */
        virtual Index_t local_cols () const
        {
            return col_end () - col_start ();
        }

        /** The index of the first matrix row stored on this processor.*/
        virtual Index_t row_start () const = 0;

        /** The index of the last matrix row (+1) stored on this processor. */
        virtual Index_t row_end () const = 0;

        /** The index of the first matrix column owned by this processor. */
        virtual Index_t col_start () const = 0;

        /** The index of the last matrix column (+1) owned by this processor.*/
        virtual Index_t col_end () const = 0;

        /** Set all entries to 0. */
        virtual void setZero () = 0;

        /** Set *this to the identity matrix. */
        virtual void setIdentity () = 0;

        /**
         * Set the element \p (i,j) to \p value.  Throws an error if the
         * entry does not exist. Zero values can be "stored" in non-existent
         * fields.
         */
        virtual void set (const Index_t& i, const Index_t& j, const TDataType& value) = 0;

        /** Fill the matrix *this with the list of triplets defined by the iterator range begin - end.*/
        template <typename InputIteratorType>
        void setFromTriplets (const InputIteratorType& begin, const InputIteratorType& end);

        /**
         * Add \p value to the element \p (i,j).  Throws an error if the
         * entry does not exist. Zero values can be "added" to non-existent
         * entries.
         */
        virtual void add (const Index_t& i, const Index_t& j, const TDataType& value) = 0;

        // /**
        //  * Add the full matrix \p dm to the SparseMatrix.  This is useful
        //  * for adding an element matrix at assembly time.
        //  */
        // virtual void add (const DenseMatrix<TDataType>& dm, const IndexMap_t& rows, const IndexMap_t& cols) = 0;

        // /**
        //  * Add \p spm to the rows and cols of this matrix (A)
        //  */
        // virtual void add (const SparseMatrix<TDataType>& /*spm*/, const IndexMap_t& /*rows*/, const IndexMap_t& /*cols*/) = 0;

        // /** The reference of matrix entry \p(i,j). */
        // virtual TDataType& operator() (const Index_t& i, const Index_t& j) = 0;

        /** A copy of matrix entry \p (i,j). */
        virtual TDataType operator() (const Index_t& i, const Index_t& j) const = 0;

        /** The number of nonzeros. */
        virtual Index_t nonZeros () const = 0;

        /**
         * Get the CSR (Compressed Sparse Row) Storage of the sparse matrix.
         * @param ia Array, size (n+1), ia(i) points to the first column index of row i in the array ja.
         * @param ja Array, size (nnz), indices of the sparse matrix A.
         * @param a Array, size (nnz), the non-zero elements of the sparse matrix
         */
        virtual void get_csr (std::vector<Index_t>& ia, std::vector<Index_t>& ja, std::vector<TDataType>& a) const = 0;

        /** Print the elements of a matrix to the given output stream. */
        virtual void print (std::ostream& out) const = 0;

        /**
         * This function creates a matrix called "submatrix" which is defined
         * by the row and column indices given in the "rows" and "cols" entries.
         * Currently this operation is only defined for the PetscMatrixBase subclasses.
         * Note: The \p rows and \p cols vectors need to be sorted;
         *       Use the nosort version below if \p rows and \p cols vectors are not sorted;
         *       The \p rows and \p cols only contain indices that are owned by this processor.
         */
        virtual void
        create_submatrix (SparseMatrix<TDataType>& submatrix, const GeneralIndexSet& rows, const GeneralIndexSet& cols) const = 0;

        /** Create a submatrix from the row/col stride index set. */
        virtual void
        create_submatrix (SparseMatrix<TDataType>& submatrix, const StrideIndexSet& rows, const StrideIndexSet& cols) const = 0;

        /**
         * Similar to the above function, this function creates a \p
         * submatrix which is defined by the indices given in the \p rows
         * and \p cols vectors.
         * Note: Both \p rows and \p cols can be unsorted;
         *       Use the above function for better efficiency if your indices are sorted;
         *       \p rows and \p cols can contain indices that are owned by other processors.
         */
        virtual void create_submatrix_nosort (SparseMatrix<TDataType>& /*submatrix*/,
                                              const GeneralIndexSet& /*rows*/,
                                              const GeneralIndexSet& /*cols*/) const = 0;

        /**
         * Create a submatrix from the row/col stride-index-set.
         * \p rows and \p cols can contain indices that are owned by other processors.
         *
         */
        virtual void create_submatrix_nosort (SparseMatrix<TDataType>& /*submatrix*/,
                                      const StrideIndexSet& /*rows*/,
                                      const StrideIndexSet& /*cols*/) const = 0;

        /** Copy the diagonal part of the matrix into \p dest. */
        virtual void get_diagonal (Vector<TDataType>& dest) const = 0;

        /** Copy the transpose of the matrix into \p dest, which may be *this. */
        virtual void get_transpose (SparseMatrix<TDataType>& dest) const = 0;

        /** Copy the inverse of this matrix into \p dest, which may be *this. */
        virtual void get_inverse (SparseMatrix<TDataType>& dest) const = 0;

        /** Scale all elements of this matrix by \p scale */
        virtual void scale (const TDataType scale) = 0;

        /** Compute Y = A*X for matrix \p X. */
        virtual void matrix_matrix_mult (SparseMatrix<TDataType>& X, SparseMatrix<TDataType>& Y) = 0;

      protected:
        // bool      mInitialized = false;    //< whether or not the matrix has been initialized
        Sparsity* mSparsity    = nullptr;  //< the data pattern across processes.
    };

    /**
     *  Mimic Eigen3 SparseMatrix implementation.
     * 
     * template<typename Scalar , int Options_, typename StorageIndex_ >
     * template<typename InputIterators >
     * void Eigen::SparseMatrix< Scalar, Options_, StorageIndex_ >::setFromTriplets	(const InputIterators& begin, const InputIterators& end)
     */
    template <typename TDataType>
    template <typename InputIteratorType>
    void SparseMatrix<TDataType>::setFromTriplets (const InputIteratorType& begin, const InputIteratorType& end)
    {
        for (auto it = begin; it != end; ++it)
        {
            this->set (it->row (), it->col (), it->value ());
        }

        this->assemble ();
    }

    // // Alias for the wrappers to avoid potential conflicts.
    // template <class TDataType>
    // using TVector = Vector<TDataType>;
    // template <class TDataType>
    // using TDenseMatrix = DenseMatrix<TDataType>;
    // template <class TDataType>
    // using TSparseMatrix = SparseMatrix<TDataType>;

    /**
     * @brief Define the interfaces to create a family of algebraic objects.
     *
     * Abstract Factory Pattern is taken to provide an interface for creating families of related 
     * or dependent objects without specifying their concrete classes.
     *
     * @see Erich Gamma. Design Patterns:elements of reusable object-oriented software. Addison Wesley, 1994.
     *
     * @author nene
     * @date February, 2025.
     */
    template <typename TDataType>
    class ALGEBRA_EXPORT Space
    {
      public:
        virtual Vector<TDataType>* createVector ()
        {
            return nullptr;
        };

        virtual DenseMatrix<TDataType>*  createDenseMatrix ()
        {
            return nullptr;
        }

        virtual SparseMatrix<TDataType>* createSparseMatrix ()
        {
            return nullptr;
        }
    };

    //@}
}  // namespace Algebra
SG_NAMESPACE_CLOSE

#endif  // Algebra_Space_h